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■ Heat capacity curves as functions of temperature for classical atomic clusters 
^ ■ bound by pairwise Lennard-Jones potentials were calculated for aggregate 
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sizes 4 < < 24 using Monte Carlo methods. J-walking (or jump-walking) 
was used to overcome convergence difficulties due to quasi-ergodicity in the 
■ solid-liquid transition region. The heat capacity curves were found to dif- 

fer markedly and nonmonotonically as functions of cluster size. Curves for 
= 4, 5 and 8 consisted of a smooth, featureless, monotonic increase through- 
out the transition region, while curves for N = 7 and 15-17 showed a distinct 
shoulder in this region; the remaining clusters had distinguishable transition 
heat capacity peaks. The size and location of these peaks exhibited "magic 
number" behavior, with the most pronounced peaks occurring for magic num- 
ber sizes of = 13, 19 and 23. This is consistent with the magic numbers 
found for many other cluster properties, but there are interesting differences 
for some of the other cluster sizes. Further insight into the transition region 
was obtained by comparing rms bond length fluctuation behavior with the 
heat capacity trends. A comparison of the heat capacities with other cluster 
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properties in the solid-liquid transition region that have been reported in the 
literature indicates partial support for the view that, for some clusters, the 
solid-liquid transition region is a coexistence region demarcated by relatively 
sharp, but separate, melting and freezing temperatures; some discrepancies, 
however, remain unresolved. 
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I. INTRODUCTION 



The investigation of small atomic and molecular clusters is an active area of current 
research, both theoretically§Hll and experimentally.EUll Their small size (typically 3 to 150 
atoms) makes clusters ideal candidates for computer simulation studies, and most theoretical 
studies to date have used Monte Carlo&lfl and molecular dynamics methods-E^ E3 Because 
of the large computational demands inherent in these simulation methods, most early the- 
oretical investigations were limited to classical rare gas clusters bound by simple pairwise 
additive Lennard- Jones potentials, but continuing improvements in computer technology 
have led to an increasing number of more interesting metal cluster simulations requiring 
more sophisticated intermolecular potentials for physically realistic representation,!! and 
have made quantum simulations based on Fourier path integral methods practical.!! 

One major motivation for the study of clusters is the insight that can be provided for 
understanding the transition from finite to bulk behavior. Many cluster properties show 
a very irregular dependence on the aggregate cluster size A^. For example, the difference 
in the binding energy between successive pairs of rare gas clusters is especially marked for 
clusters of size = 13, 19, 23, 26, 29, 39, 46, 49, and 55, and dramatically smaller for their 
immediate neighbors.!! The mass spectrum intensity of Xe clusters formed by condensation 
in supersonic jets was found to increase slightly for sizes N = 13, 19, 25, 55, 71, 87, and 
147 and then decrease sharply for the clusters immediately following in size.lli Calculations 
of the Gibbs free energy of formation for argon clusters with sizes ranging from 2 to 20 
indicated local minima for = 7, 13 and 19.0 This nonuniform dependence of certain 
cluster properties on size, often termed "magic number" effects, has been of great theoretical 
interest. Much work has been done in relating magic number sequences to cluster structure 
in terms of "soft" sphere packings, since many magic number sizes such as iV = 13, 19, 55, 
and 147 correspond to compact icosahedral based structures.ii^0 

An especially important phenomenon exhibiting magic number effects is cluster "melt- 
ing." Because of their small size, clusters do not have the sharp first order transition charac- 
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teristic of bulk melting, but instead show transitions occurring over a range of temperaturesE 
The transition range generally decreases as the cluster size increases, coalescing to the bulk 
transition temperature as approaches infinity, but it shows a wide variety of behavior 
for smaller sizes. Berry and coworkersS have postulated that for certain sized clusters, the 
transition region is in fact a coexistence region where both "solid-like" and "liquid-like" 
isomers dynamically coexist, fiuctuating back and forth between relatively long lived "solid" 
and "liquid" states. This implies that these clusters have sharp, but different, "freezing" 
and "melting" temperatures. Some of their molecular dynamics studies revealed that clus- 
ters with sizes = 7, 9, 11, 13, 15, and 19 have especially high melting temperatures and 
pronounced coexistence ranges, while clusters with sizes = 8, 14 and 20 melt at much 
lower temperatures and exhibit no coexistence.^ 

The problem of cluster phase transitions, has been one of great interest, and one not lack- 
ing in controversy (Ref. ^ provides a good review). While many studies are consistent with 
the view of a dynamic coexistence between solid-like and liquid-like isomers in the transition 
region, other studies emphasize, instead, a relatively smooth progression of isomerizations 
between locally similar but globally distinct configurations occurring throughout the tran- 
sition region.0 In addition, it is becoming increasingly clear that the problem of incomplete 
sampling of configuration space for simulations run in the transition region is much more 
serious than had been appreciated earlier, and that extremely long simulations (much longer 
than those routinely used even a few years ago) are necessary to prevent large systematic 
errors in certain properties that can result in misleading interpretations. Given the large 
amount of contradictory evidence that has been accumulated over the years, as well as the 
problem of not knowing how reliable some of it is, it is apparent that the nature of cluster 
melting remains still very much unresolved. 

There is no problem concerning cluster behavior at extreme temperatures. At very low 
temperatures, clusters are locked into solid structures having near rigid rotations and small 
amplitude, harmonic-like vibrations; there is no conversion between different isomers. At 
sufficiently high temperatures, the clusters dissociate, implying that for simulations run 



under free volume conditions, the average energy vanishes in the limit of an infinite number 
of configurations.EUll Consequently, it is usual in Monte Carlo simulations to confine the 
cluster within a perfectly reflecting constraining potential having a typical radius of a few 
atomic radii that is centered on the cluster's center of mass.0 Under such conditions, the 
high temperature limit corresponds then to a highly compressed vapor confined to a small 
spherical cavity. This limiting behavior is evident in Fig. ||, which shows the average internal 
energy and heat capacity as functions of temperature for a typical small classical rare gas 
cluster. For T — 0, the energy approaches the global minimum potential energy, and 
the heat capacity approaches the classical equipartition of energy result for a harmonically 
bound nonlinear system, 3A^ — 3. For T — oo, the internal energy and heat capacity both 
approach their ideal gas limits. 

What happens in between is less clear. In this case, the heat capacity curve has two 
distinct peaks, which would indicate two "phase" changes. Temperature driven first-order 
transitions in small systems are typically spread out over a finite range of temperatures, 
with the size of the transition region corresponding to the width of the heat capacity peak, 
and the transition temperature corresponding to the heat capacity maximum.!! The larger, 
higher temperature heat capacity peak can be therefore associated with the transition from 
a condensed phase to the vapor phase (that is, cluster dissociation), while the smaller, 
lower temperature peak represents a solid-liquid transition. Not surprisingly, the size and 
location of the higher temperature peak have a sensitive dependence on the cluster simulation 
boundary conditions, and thus vary considerably as functions of the constraining radius. 
Their dependence on the cluster aggregate size is less pronounced, though, and quite regular. 
Hence the information conveyed in this region is somewhat artificial and lacking in physical 
importance. This is not the case for the solid-liquid transition, however. This peak is 
largely independent of the constraining radius, provided the radius is not so small that the 
free rearrangement of the cluster is limited. In addition, the solid-liquid transition region 
has a very marked dependence on cluster size; while all clusters confined in sufficiently large 
constraining spheres show a large dissociation peak in their heat capacity-temperature curve, 
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not all clusters show peaks in the solid-liquid transition region, and for those that do, the size 
and shape of the peak have a very nonuniform dependence on cluster size. Does the presence 
of a heat capacity peak in the solid-liquid transition region then necessarily imply a solid- 
liquid transition temperature corresponding to the peak temperature, and concomitantly, 
does the absence of a heat capacity peak there rule out a definite solid-liquid transition 
temperature? What is the connection between the transition region as defined by the heat 
capacity peak width and the coexistence region defined by the sharp, yet distinct, freezing 
and melting temperatures obtained from some molecular dynamics simulations? One would 
expect the two ranges to mostly coincide, but as will be seen, this is not always the case. 

Despite the importance of heat capacities in characterizing cluster phase transitions, few 
cluster heat capacity studies have been reported in the literature, and most of these have 
been limited to special cases, particularly = 130'&B and = 55.00 This is due to the 
notorious difficulty in extracting heat capacities with sufficient accuracy from simulations. 
A typical example of the poor convergence of heat capacities is given in Fig. [1|. The heat 
capacity curve obtained from averages taken over 10^ Monte Carlo passes (or 23 x 10^ moves) 
for each temperature point (a simulation length that was typical ten to fifteen years ago) is 
not even qualitatively correct, with the solid-liquid transition peak completely absent. The 
curve obtained from walk lengths of 10^ passes (a length routinely achieved five years ago) 
does show a small solid-liquid transition peak, but one that is very much diminished. Even 
the curve based on 10^ passes is quantitatively poor in this region. This poor convergence 
behavior is a result of the incomplete sampling of configuration space by the random walker, a 
phenomenon denoted as "quasi-ergodicity" by Valleau and Whittington.0 Quasi-ergodicity 
arises in systems where the sample space contains two or more regions that have a very 
low transition probability, resulting in bottlenecks that effectively confine the sampling to 
only one of the regions.0 In clusters such as Ar23, the potential hypersurface consists of 
very deep wells, corresponding to stable solid isomers, separated by very large barriers. 
The resultant dichotomy of time scales characterizing the Monte Carlo random walks (or 
likewise, the molecular dynamics time evolution) produces rapid small-scale motion within 
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the wells, but very slow large-scale movement between them.El Because of the finite length 
of the walks (or trajectories), systematic errors result that diminish only with increasing 
length, disappearing in the limit of infinitely long walks (or trajectories). Interestingly, it is 
precisely these conditions that Berry and coworkers claim are necessary for the existence of a 
solid-liquid coexistence region in clusters;^ it appears that the very nature of this region that 
makes it so physically intriguing also conspires to make its elucidation exceedingly difficult. 

Recently, new methods have been developed that substantially reduce quasi-ergodicity 
and dramatically increase convergence in cluster Monte Carlo simulations in the transition 
regions, thus allowing even computationally difficult properties such as cluster heat capacities 
to be obtained to high accuracy. These methods are characterized by their emphasis on large 
scale motions, or jumps, of the random walker to other regions of configuration space and 
thus have been collectively termed "J- walking" (for jump- walking) .13 J- walking has been 
successfully applied to several systems, including binary metal clusters,!! water clusters,!! 
and clusters absorbed on surfaces.0 Tsai and JordaniiS have combined J-walking with 
histogram methods,^ and recently, J-walking was extended to quantum systems. @ With 
J-walking, it is now feasible to systematically study the dependence of the heat capacity on 
cluster size for a fairly large range of cluster sizes; the solid curves in Fig. |T] were obtained 
using J-walking. 

I begin in Section |I| with a brief review of the J-walking method and a summary of 



the calculations undertaken. Section |ITl| describes the results obtained using J-walking for 
the heat capacities and potential energies as functions of temperature and cluster size for 
Lennard- Jones clusters ranging in size from 4 to 24 atoms. Also included are results obtained 
from similar standard Monte Carlo simulations, including rms bond length fiucuations, as 
well as comparisons with molecular dynamics simulations and other calculations reported 
in the literature for rare gas clusters in this range. Finally, in Section 0, I discuss my 
findings in the context of some of the controversies and ambiguities concerning cluster phase 
transitions that remain unresolved. 
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II. METHOD 



J-walking addresses the problem of quasi-ergodicity in standard Monte Carlo simula- 
tions based on the sampling algorithm proposed by Metropolis et al.^ hereinafter called 
Metropolis Monte Carlo, by coupling the usual small-scale Metropolis moves with occasional 
large-scale jumps that move the random walker to different regions of configuration space. 
For classical systems governed by a potential V^(r), the Metropolis method uses a random 
walk to sample configuration space, with each attempted move from the current location 
to a trial location rj determined by the acceptance probability, 

P = min{l,g(r/|ri)}, (1) 

where 

^^^^"^^^ - T(r,|r.)p(r.)' 

is the acceptance ratio, p(r) = exp{—f3V{r)} is the Boltzmann distribution with Z the 
standard configuration integral and (3 = 1/A;bT the temperature parameter, and T(r'|r) is 
the transition matrix or sampling distribution. The sampling distribution is usually gener- 
ated from uniform deviates ^ over a finite stepsize range A to giveii 

T( I ^ 1^/^ |r,-r,|<A/2, 

T{rf\ri) = ( (3) 
I otherwise, 

and thus 

g(r;|r,) = exp{-/3[y(r^) - ^(r,)]}. (4) 

Attempted moves are generated according to rj = rj + (,^ — |)A, with the maximum stepsize 
A/2 usually adjusted to give acceptance probabilities of approximately 50%. The required 
size depends on the width of the Boltzmann distribution, and thus decreases with increasing 
p. This temperature dependence can lead to quasi-ergodic behavior whenever the stepsize 
becomes too small relative to the potential barrier heights and widths; the walker becomes 
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effectively trapped within a region of configuration space if too many small steps are required 
to surmount large, steep barriers within the duration of the walk. Thus, for a given walk 
length, there is a threshold temperature such that walks undertaken below that temperature 
are quasi-ergodic.0 

The Boltzmann distribution's dependence on temperature that can lead to quasi- 
ergodicity can also be exploited to overcome it. Although the form of the Boltzmann 
distribution is largely dependent on the form of the underlying potential, with distribu- 
tion maxima corresponding to potential minima, the widths of the distribution peaks have 
an inverse dependence on /5, with higher temperatures resulting in wider distributions. In 
essence, higher temperature walkers are less constrained because their larger stepsizes allow 
them to overcome the potential barriers more effectively. Hence, a Boltzmann distribution 
generated at some higher temperature where the sampling is ergodic can be used by a lower 
temperature walker to make large-scale moves in configuration space. J-walking then occa- 
sionally replaces the usual attempted Metropolis moves with attempted jumps to positions 
occupied by such a higher temperature walker (J- walker). Because the peaks in the Boltz- 
mann distribution correspond to the potential minima, the J-walker's motion remains biased 
about the minima, greatly increasing the likelihood an attempted jump would be accepted. 

This scheme is equivalent to replacing the usual Metropolis transition matrix given by 
Eq. d^) with the Boltzmann distribution at the higher temperature whenever a jump is 
attempted, 

Tj(rj|r,) = Z-^ exp{-/?jV^(rj)} for < O < Pj. (5) 

where Pj is the jump attempt probability, fij is the J-walker temperature parameter, and 
rj is the J-walker's location. This gives an acceptance ratio of 

qj{jj\r>) = exp{(/5j - P)[V{vj) - \/(r,)]}. (6) 

The likelihood of an attempted jump being accepted depends on the two Boltzmann distri- 
butions' overlap, which is reflected by the value of qj. In the high temperature limit (3j — >■ 0, 
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the acceptance ratio reduces to the standard Metropohs expression given in Eq. (H). Because 
the Boltzmann distribution broadens as the temperature increases, J-walking in this hmit 
reduces to indiscriminate jumping with a very large stepsize, and so the probabihty of a 
jump being accepted becomes very small. In the limit (3j — > /5, qj{vj\Yi) 1 since the low 
temperature walker is now effectively sampling from its own distribution. 

Two complementary implementations for generating the classical J-walker Boltzmann 
distributions were originally presented.^ The two are another example of the commonly 
encountered trade-off between computer memory and speed. The first ran the J-walker 
in tandem with the low temperature walker, with the low temperature walker occasionally 
attempting jumps to the current J-walker position simply by using the current J-walker coor- 
dinates as its trial position. In the second implementation, the J-walker was run beforehand 
and the configurations generated during the walk stored periodically in an external array. 
Subsequent jump attempts were made by accessing the stored configurations via randomly 
generated indices. The tandem walker scheme has very little memory overhead, but has the 
disadvantage of requiring that the J-walker be moved an extra number of steps whenever a 
jump is attempted in order to reduce correlations between the two walkers that otherwise 
result in systematic errors. The number of extra steps needed depends on the temperature 
difference between the two walkers, increasing the computation time greatly as the differ- 
ence becomes larger. This implementation then is more suited for parallel computers. On 
the other hand, the use of external distributions has only a modest computational overhead 
(mostly the time required to generate the distributions), but requires very large storage fa- 
cilities for handling the distribution arrays. Because fast workstations having tens of Mb of 
RAM and Gb of disk storage are now affordable and quite common, while access to parallel 
computers is still much more limited, the implementation using externally stored distribu- 
tions remains the method of choice and was the one used for all the J-walking calculations 
reported in this study. 

J-walking simulations were run for all clusters in the range 4 < < 24. The clusters 
studied were modeled by the usual pairwise additive Lennard- Jones potential, 
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Small clusters are known to become unstable beyond a threshold temperature that varies 
with the cluster size.0'0 For the Lennard- Jones potential under free volume conditions, the 
average energy vanishes in the limit of infinitely long walks. Consequently, the choice of 
boundary conditions can greatly affect some cluster properties at higher temperatures.^ I 
have followed Lee, Barker and Abraham^ and have confined the clusters by a perfectly 
reflecting constraining potential of radius Rc centered on the cluster's center of mass. To 
maintain a common set of boundary conditions throughout the survey, the constraining 
radius was identical for all the clusters studied, with = ^(J. 

For physically realistic systems such as clusters, the Boltzmann distributions are too 
narrow for a single distribution to be used for J-walking over the entire temperature do- 
main from the solid region to the dissociation region, and so the distributions had to be 
generated in stages. For each cluster, an initial J-walker distribution was generated from a 
long Metropolis walk at a temperature high enough for the sampling to be ergodic. This 
distribution was then used for J-walking runs to obtain averages of the potential energy and 
the heat capacity for a series of lower temperatures. The width of the J-walker distribution 
limited the effective temperature range for subsequent J-walking simulations. When the 
temperature difference between the distribution and the low temperature walker became 
too large, very few attempted jumps were accepted since the low J-walker configurations 
likely to be accepted were in the low energy tail of the distribution. At the temperature 
where the jump acceptance started becoming too small,^ a new distribution was generated 
from the previous one, using J-walking to ensure it was ergodic as well. This distribution 
was then used to obtain data for the next lower temperature range, and then to generate 
the next lower temperature J-walker distribution, and so on, until the entire temperature 
domain was spanned. To ensure that the higher temperature J-walker distributions were 
free of any systematic errors large enough to adversely affect the sampling, the temperature 
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ranges were occasionally overlapped so that J-walking data collected near the end of one 
distribution's useful temperature range (as indicated by the low jump acceptance) could 
be compared with the data obtained from the beginning of the next distribution's useful 
temperature range (where the jump acceptance rate was high). For example, a J- walker 
distribution generated at T = 15 K could be used to generate data over the temperature 
range 11<T<15K. By generating the next distribution at T = 12 K for use over the 
range 8 < T < 12 K, data for the overlapping range 11 < T < 12 K obtained from both 
distributions could be compared and checked for consistency. All the distributions checked 
in this manner were found to be consistent. 

The classical internal energy and heat capacity were calculated by the usual expressions 
for an A^-atom cluster (in reduced units), 

(U*) = ^- + {V*) , (8) 
3N {{V*)^)-{V*f 

\^V) - ^ + ^J^2 ' W 

where U* = U/e,V* = V/e, Cy = Cv/ks, and T* = ksT /e. In order to more easily compare 
the classical results with quantum simulations to be undertaken later, curves for (f/*) and 
{Cy) as functions of T* were actually generated using Ne parameters, with a = 2.749 A and 
e = 35.6 K. The temperature mesh size was typically AT = 0.02 K, and was reduced to 
0.01 K for the solid-liquid transition regions. Similarly, J-walker distributions were generated 
at Ne temperatures. All the J-walking simulations run for each temperature consisted of 
10^ warmup passes followed by 10"^ passes with data accumulation (100 walks of 10^ passes 
each); the jump attempt frequency was Pj = 0.1 throughout. 

Table | lists the J-walker distribution particulars for clusters in the range 4 < A^ < 12, 
while Table ^ lists them for the range 13 < A^ < 24. For each cluster, the initial J-walker 
distributions were obtained at temperatures corresponding to the high temperature side of 
the cluster dissociation peak, where all the cluster configurations were completely fluid. In 
the original J-walking study for A^ = 13,0 the initial J-walker distribution was obtained at a 
temperature corresponding to the liquid region (just past the solid-liquid transition region) 
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where the samphng was thought to be ergodic, but in fact turned out not to be fully ergodic. 
When J-walking was extended to quantum systems and tested on = 7 clusters,El it was 
discovered that there was also much quasi-ergodicity in the region corresponding to the low 
temperature side of dissociation peak, and so it was necessary to begin the temperature 
scans at an even higher temperature. This region corresponds to a coexistence region where 
both liquid-like and partially dissociated configurations are in dynamic coexistence, and 
sampling difficulties can arise when single atoms leave the cluster and wander about nearby 
for several passes before recombining with the cluster. Effects of quasi-ergodicity in this 
region can be seen in the heat capacity curves for Ar23 in Fig. |l[ On the high temperature 
side of the dissociation peak (T ^ 80 K) the J-walker curve coincides with the Metropolis 
curves for both 10^ and 10^ passes, while on the low temperature side of the peak (50 K 
^ T < 75 K ) there are substantial discrepancies between the J-walking curve and the 
Metropolis curve obtained from 10^ passes, and smaller discrepancies between the J-walking 
curve and the Metropolis curve obtained from 10^ passes. Although these discrepancies are 
much smaller than those appearing in the solid-liquid transition region, it is important to 
minimize systematic errors in the J-walker distribution since the sensitivity of these errors 
on subsequent J-walking runs is not known. 

All J-walker distributions consisted of 10^ configurations, except for A^ = 4 and A^ = 8, 
which had 1.25x 10^ and 0.75x 10^ configurations, respectively. The configuration energy and 
center of mass coordinates were also stored with the cluster coordinates so that they would 
not have to be recalculated during subsequent J-walker simulations whenever a jump was 
accepted. The required storage space for these distributions varied from 122 Mb (for A^ = 4) 
to 580 Mb (for A^ = 24), and so they were stored in 10 or 50 separate files, depending on the 
amount of computer primary memory that was available.0 During a J-walking run, a file was 
randomly selected and read into memory. This array was used for sampling jump attempts 
for a while before being randomly replaced by another file. By making the files as large 
as possible within the memory constraints, disk access was reduced during the J-walking 
simulations, greatly improving their efficiency. Because the configurations generated during 
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Metropolis walks are highly correlated, they were stored in a periodic fashion. Initially, 
almost all distributions were generated by storing configurations every 100 passes, as can 
be seen in Table |I| for the smaller clusters. As the cluster size was increased, generating 
such distributions became very time consuming. Not only did it take longer to generate 
a distribution for a larger cluster, but more distributions were required to span the entire 
temperature range. Thus, for the larger clusters, only the widest distributions containing the 
greatest variety of configurations were sampled every 100 passes. Narrower distributions were 
sampled every 50 passes, and very narrow distributions for low temperatures corresponding 
to the solid region were sampled every 10 passes. Fig. |^ shows potential energy histograms 
for = 11 and = 19 clusters. The smaller cluster has a potential energy range about 
half that of the larger cluster, and so only five J-walker distributions were needed to span 
the entire temperature range from the solid region to dissociation, compared to the ten 
distributions required for the larger cluster. The spacing of the distributions requires some 
judgment. On the one hand, they should be spaced as far apart as possible to minimize the 
computational time needed to generate them. On the other hand, they should not be spaced 
so far apart that the number jumps accepted becomes too low and systematic errors arising 
from quasi-ergodicity corrupt the distributions. As can be seen in the plots, substantial 
overlap between adjacent distributions was required for reasonable jump acceptances (these 
are also listed in Tables | and ||) . 

Correlations in the J-walker distributions were further reduced by writing the distribution 
files in a parallel fashion. All output distribution files were opened at the start of the program 
and configurations were written to each in turn, rather than writing to each file sequentially, 
one at a time. Thus, each distribution file contained configurations sampled from the entire 
run. Unfortunately, this resulted in highly fragmented files, which greatly increased the 
time required to later read a distribution into memory, since the files were all physically 
interleaved on the disk drive. This problem was easily solved, however, by simply copying 
the distribution files to a different directory, where they were written to the disk sequentially. 

Since J-walking is still new, and experience with it is still being obtained, standard 

14 



Metropolis simulations were also run for each cluster. These provided a check of the J- 
walking results for those temperature regions where quasi-ergodicity in the Metropolis runs 
was not a problem, as well as revealing trends in the systematic errors arising from quasi- 
ergodicity in the Metropolis runs as functions of cluster size. Temperature scans were also 
generated using Ne parameters, with the temperature mesh size AT = 0.02 K. For each 
temperature, simulations consisted of 10^ warmup passes, followed by 10'' passes with data 
accumulation (again, 100 walks of 10^ passes each). The scans were started at T = 0.2 
K from the global minimum configuration,iii2l and continued past the cluster dissociation 
peak. The final configuration for each temperature was used as the initial configuration for 
the next temperature. For all of the clusters examined, the Metropolis and J-walking results 
agreed qualitatively throughout the entire temperature ranges, and agreed quantitatively 
throughout most of the temperature ranges, with discrepancies occurring mostly in the 
transition regions affecting the peak heights, especially for the larger sizes; the discrepancies 
seen in Fig. |1| for Ar23 were the largest obtained in this study. 

One of the few cluster properties not amenable to calculation using J-walking is the 
relative rms bond length fluctuation. 



The averages are taken over an entire walk, and so are dependent on the walk length. With 
J-walking, a "walk" is really a collection of short excursions of a few steps duration beginning 
from the different points in configuration space that are accessed whenever a jump attempt 
is accepted, and so there is no continuous evolution of i-j bonds over the entire length of 
the walk. 

Rms bond length fluctuations have been calculated in many previous molecular dynam- 
ics and Metropolis Monte Carlo simulations and have been used as indicators of solid-liquid 
phase changes.0""i0i3 The fluctuations typically undergo a gradual, nearly linear increase as 
the temperature is increased in the solid region, then rise sharply in the transition region as 
the clusters begin to acquire enough energy to overcome the potential barriers to rearrange- 




(10) 
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ment, and then level off somewhat in the liquid region. For magic number clusters such as 
= 13, the rise can be very abrupt, while for other clusters, it is more gradual.ll Like the 
heat capacity, 6{T) is very difficult to calculate accurately in the transition region because 
of its sensitivity to quasi-ergodicity. This is evident in Fig. ^, which shows averages of S{T) 
for = 13 obtained from 10 Metropolis walks of lO'', 10^ and 10^ passes in length; the 
lack of convergence in the transition region is clear. The systematic error throughout the 
transition region results in the S{T) curves for 10^ and 10^ passes occurring too high in the 
temperature region spanned by the heat capacity peak, whereas the curve for 10^ passes rises 
sharply just as the heat capacity curve begins its rise. The curve for 10^ passes is similar to 



the curve reported in Ref. which was obtained from Metropolis walks of a similar length. 
This behavior can also be understood in terms of the potential energy hypersurface govern- 
ing the random walks. For = 13, the deep potential wells associated with various stable 
isomers are separated by large barriers. Temperatures near the point T/ where S{T) rises 
sharply correspond to energy distributions having a small fraction of energies in their tails 
that are on the threshold of overcoming the barriers to isomerization. Thus a walk length 
of 10^ passes that resulted in even a few barrier crossings would be averaging the different 
fluctuations for different isomers, resulting in much larger values than those obtained from 
walks of 10^ or 10^ passes where the walker remained trapped in a single isomeric form for 
the duration of the walk. 

This great sensitivity of the S{T) curve on the walk length brings into question then the 
appropriateness of applying to clusters the Lindemann criterion of liquid behavior occurring 
for 6 > 0.1, which has its basis from the study of bulk systems. For clusters, a value of 
(5 ~ 0.1 obtained from sufficiently long walks actually corresponds to the threshold of iso- 
merization where the cluster spends most of its time in some solid-like form with occasional 
isomerizations to another solid-like form, rather than to liquid behavior where the cluster 
atoms readily undergo diffusion. The region following the sharp rise where the curve levels 
off is more representative of liquid behavior in clusters. This point has been previously men- 
tioned in accounting for the different energies observed in molecular dynamics simulations 
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between the onset of large bond length fluctuations and the onset of diffusion obtained from 
power spectra.il is also much lower than the value obtained from molecular dynamics 
simulations reported in Ref. 0, although that value may also be too high considering the 
great sensitivity of 6 on the simulation length found in the Monte Carlo simulations. Care 
must be taken when comparing S values obtained from molecular dynamics simulations with 
those obtained from Monte Carlo simulations, though. Average rms bond length fluctua- 
tions can depend greatly on the ensemble used, and discrepancies between results obtained 
for Ari3 in the microcanonical ensemble with molecular dynamics and those obtained in the 
canonical ensemble with Monte Carlo have been noted previously.ii Despite these problems, 
S{T) curves can still be very useful for providing insight into the solid-liquid transition region 
since their behavior also has a sensitive dependence on the cluster size. 

Rms bond length fluctuations were not originally calculated with the Metropolis sim- 
ulations described above, but in order to obtain more insight into some discrepancies in 
the transition region between the J-walking results and results obtained from molecular dy- 
namics calculations reported in Ref. ^ that became apparent after the original Metropolis 
simulations had been run, temperature scans of 6 were subsequently obtained from addi- 
tional Metropolis simulations. These were calculated at each temperature as averages taken 
over 10 walks, each having a length of 10^ passes. 



III. RESULTS 

Curves of the heat capacity per atom as functions of temperature for clusters ranging 
in size from = 4 to 12 and 13 to 24 are shown if Figs. ^ and |^, respectively. These 
were obtained from the J-walking studies. As is evident in these figures, the curves differ 
markedly and nonmonotonically with cluster size. There are three kinds of curves. For 
clusters of size A^ = 4, 5 and 8, the heat capacity rises slowly from its classical equipartition 
value at T* = until the temperature approaches the cluster dissociation region, where it 
then rises sharply. There is no peak in the solid-liquid region, although a very weak shoulder 
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is barely perceptible. Curves for clusters of size N = 7, and 15-17 also have no peak in 
the solid-liquid transition region, but they do show a distinct shoulder characterized by two 
inflection points, one at a lower temperature T£, about halfway up the shoulder, and the 
other at a higher temperature T^, at the top of the shoulder. All the other clusters studied 
have discernible heat capacity peaks in this region that can be characterized by an inflection 
point on the low temperature side of the peak (at T£), the peak maximum itself (at Tp), 
and a local minimum on the high temperature side of the peak (at T^) between the peak 
and the dissociation region peak. These temperatures are indicated in Figs. § and ^ by 
the large ticks on the temperature axes; their values, together with their corresponding heat 



capacities, are also listed in Table |ITT[ The values for the peak and local minimum in each case 
were obtained by smoothing the J-walker data (note, however, the curves shown in Figs. ^ 
and ^ are the raw, unsmoothed J-walker data), interpolating the smoothed data to obtain 
a finer mesh size, and then searching for the maximum and minimum, respectively. The 
inflection points were obtained by using finite differences in the smoothed, interpolated data 
to generate derivatives, and then searching the resulting derivative curves for the maximum 
and minimum. Because there was some variation in the values obtained from various fits 
generated with different smoothing parameters, several fits were done for each cluster; the 
values reported are the averages obtained from the fits, and the uncertainty estimates are 
the standard deviations.0 

For those clusters having peaks in the solid-liquid region, the transition temperature can 
be associated with the peak maximum, and the transition range with the peak width,S while 
for those clusters having shoulders in this region, the transition temperature can be estimated 
as the inflection point at the top of the shoulder. Although the transition temperature 
thus defined can be obtained easily for each curve, defining the transition range is more 
ambiguous. While the inflection points and local minima that characterize each peak may 
be unambiguous, their relation to the width of the transition region is less clear. In their 
original J-walker study of classical systems, Frantz, Freeman and Dollil observed that the 
onset of the transition region in Aria could be characterized by a sudden, sharp rise in the 



heat capacity standard deviation obtained from Metropolis Monte Carlo simulations. For 
very low temperatures corresponding to the Ari3 solid region, the standard deviation was 
very low, and increased slowly with increasing temperature until about T* = 0.20, where 
it quickly rose about tenfold, reaching a maximum at T* = 0.25, and then decreased to a 
minimum at T* = 0.34 before rising again in the dissociation region. If the temperature 
of the sharp rise in the heat capacity standard deviation, T*, was taken to signify the 
beginning of the transition region, and the latter temperature was taken to signify the end 
of the transition region, then the resulting temperature range was found to coincide very well 
with the Ari3 transition region obtained from molecular dynamics calculations by Berry and 
coworkers.0 Moreover, this behavior of the heat capacity standard deviations was consistent 
with their interpretation of the transition region being a coexistence region marked by sharp, 
but unequal, melting and freezing temperatures, and Tj; below Tj, only solid-like isomers 
exist, and above T^, only liquid-like forms exist, but in between in the coexistence region, 
both forms dynamically coexist. Thus, the fluctuations in the heat capacity were very small 
for low temperatures since the clusters were locked into their minimum energy configurations, 
undergoing small, mostly harmonic oscillations. At the freezing temperature, the motion 
was highly anharmonic as the clusters began to have sufficient energy to overcome the large 
barriers in configuration space separating the deep wells associated with the different isomers. 
This hindered isomerization led to the sharp increase in the heat capacity fluctuations since 
the random walks were mostly confined to individual wells, but occasionally escaped to other 
wells. The fiuctuations decreased again at the melting temperature because the energy was 
high enough for unhindered isomerization to take place; the random walker had free access 
to all of configuration space. To better see the connection between the solid-liquid transition 
region and the heat capacity fiuctuations, I have also included in Figs. ^ and ^ standard 
deviation curves for the heat capacities obtained from the Metropolis Monte Carlo walks. 
These curves are too noisy to quantitatively define a solid-liquid transition range, but their 
general features provide insight into this region. The arrows in Figs. ^ and ^ indicate the 
temperatures T*. By taking T* to signify the beginning of the transition region and to 
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mark the end, a definition of tlie transition range based solely on heat capacity behavior can 
be made, 

AT* = - T;. (11) 

The dependence on cluster size for T*, as well as for the heat capacity peak characteristic 
temperatures, T£, Tp and T^, is illustrated in Fig. Also shown is the dependence on 
size for their corresponding heat capacities. Because much of the behavior of cluster heat 
capacities is influenced by the structures of the lowest energy isomers. Fig. ^ also shows the 
size dependence on the differences in the minimum potential energy configurations between 
successive clusters, or the binding energy differences,^ 

AEt{N) = -AK.in(iV) = -[Knin(iV) " Knin(iV - 1)]. (12) 

The binding energy differences show magic number behavior for N = 7, 13, 19 and 23, as 
evidenced by these sizes being local maxima. 

Berry and coworkers have also investigated many of the clusters included in this study. 
In one particular molecular dynamics study of coexistence behavior in rare gas clusters, they 
reported results for clusters of size N = 7, 8, 9, 11, 13, 14, 15, 17, 19, 20, 22, 26 and 33.1 
By checking for the onset and disappearance of a bimodal distribution in the short-time 
averaged kinetic energies, they found coexistence ranges for = 7, 9, 11, 13, 15 and 19. 
They also reported temperatures md corresponding to long-time averaged kinetic energy 
values at which the rms bond length fluctuations 5 rose sharply. For comparison, I have 
superimposed these values in Figs. ^ and |, as well as the coexistence ranges AT^* = — Tj' 
for those clusters where the authors found coexistence behavior. An interesting feature of 
the molecular dynamics rms bond length fluctuations for those clusters showing coexistence 
behavior is evident in Figs. ^ and |^. For the magic number sizes A^ = 7, 13 and 19, T/ 
was near the middle of the coexistence range, nearly as high as the peak temperatures Tp, 
while for the other clusters showing coexistence behavior, A^ = 9, 11, and 15, it was near 
the beginning of the range. 
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A. Pre-icosahedral: 4 < TV < 12 



The results of the three studies (the J-walking heat capacity curves, the Metropohs Monte 
Carlo heat capacity fluctuation curves, and the molecular dynamics coexistence ranges) are 
mostly consistent for the smaller clusters with sizes less than = 13. For = 4 and 
5, there is no abrupt change in either the heat capacity or the heat capacity fluctuations 
until the beginning of dissociation region, where both rise sharply. = 6 is the smallest 
cluster to show a peak (albeit a small one) in the heat capacity curve, and the heat capacity 
fluctuations show the sharp rise and slow decrease characteristic of a solid-liquid coexistence 
region. = 7 is the smallest cluster having pentagonal symmetry in its lowest energy 
isomer ,0 and is the smallest of the magic number sizes (although it is a rather "weak" 
magic number in comparison to A^ = 13, 19 and 23), and so one might expect a more 
pronounced heat capacity peak than found for A^ = 6. However, N = 7 has only a shoulder, 
although the heat capacity fluctuations are similar to those for A^ = 6, implying a similar 
solid-liquid coexistence region. Ref. ^ also reported a coexistence range for N = 7, and as 
can be seen in Fig. ^ this range coincides quite well with the heat capacity transition range 
deflned by Eq. |ll]. The results for A^ = 8 are very similar to those obtained for A^ = 4 
and 5, and are consistent with those reported in Ref. § Wales and Berry! argued that the 
striking difference in coexistence behavior between N = 7 and A^ = 8 can be explained in 
terms of the differences in their potential energy surfaces. For N = 7, the global minimum 
corresponds to a very compact, stable structure, the pentagonal bipyramid, which lies much 
lower in energy than the other low-lying isomers and is separated from them by large barriers, 
preventing any easy isomerization. For A^ = 8, however, the global minimum corresponds 
to a pentagonal bipyramid with the eighth atom located on one of its ten equivalent faces. 
Degenerate isomers corresponding to moving the eighth atom to another face are separated 
by low barriers, as is a second isomer, which lies only slightly higher in energy. Thus N = 8 
clusters have easy access to the various potential wells and exhibit facile rearrangements. 
Heat capacity fluctuations remain small until the dissociation region since the random walker 
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experiences no bottlenecks in configuration space. Likewise for = 4 and 5. Altliougli 
tfieir minimum energy configurations (tetrahedron and trigonal bipyramid, respectively) are 
compact and stable, their barriers to interconversion to other isomers are small. 

The heat capacity results obtained for = 8 appear to be contrary to some of the results 
obtained by Adams and Stratt,@ who performed an instantaneous normal mode analysis on 
Ary, Arg and Aria. They observed transition regions for both Ary and Arg based on the 
temperature dependence of the average percentage of imaginary frequencies, but found no 
evidence of a coexistence region in either. They also reported other dynamical similarities 
between the two clusters, and interpreted their results to imply that the solid-liquid tran- 
sition should not be viewed primarily in terms of jumps between rigid and fluid structures, 
but as being a smooth progression through an ever increasing number of structural isomers 
that are locally similar, but globally distinct. The striking difference between the heat ca- 
pacity curves for A^ = 7 and 8, as well as between their heat capacity fluctuations, seems 
to indicate an inherently different transition behavior for these two clusters, and the consis- 
tency between these results and the molecular dynamics results lends further support to the 
coexistence view. However, not all the properties show such a large degree of dissimilarity. 
Self diffusion constants calculated by Adams and Stratt for Aii and Arg as functions of tem- 
perature were very similar,^ and curves of the rms bond length fluctuations 5 as functions 
of temperature for these two clusters are more alike than different. Although the featureless 
nature of the N = 8 heat capacity curve might imply that this cluster does not have a 
solid-liquid transition region, the diffusion constant behavior and the rms bond length fluc- 
tuation behavior suggest otherwise. The ^(T) behavior can be seen in Fig. |^, which shows 
curves for A^ = 4-6 and 7-9. Both the N = 7 and 8 curves show the sharp rise in S{T) 
that is characteristic of a transition. Although the heat capacity curves for A^ = 4, 5 and 8 
are very similar in their lack of significant features in the solid-liquid transition region, only 
the A^ = 4 S{T) curve has a likewise featureless rise through the transition region, implying 
that no solid-liquid transition occurs before the cluster dissociates. Surprisingly then, the 
absence of a peak or shoulder in the heat capacity curve does not necessarily correspond to 
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the absence of a transition. One important difference for the N = 7 S{T) curve is that its 
sharp rise occurs at a significantly higher temperature than that of its immediate neighbors 
(it is almost as high as that for = 9). This was also the case for the molecular dynamics 
results reported in Refs. ^ and ^ As will be seen, this is a common feature of the other 
magic number clusters, and is another manifestation of their inherent stability. 

The heat capacity curves for = 9 to 12 all exhibit small heat capacity peaks in the 
solid-liquid transition region, and their Metropolis Monte Carlo heat capacity fluctuations 
show the sharp rise and slow decrease characteristic of coexistence behavior. Ref. ^ reported 
results for A^ = 9 and 11, and coexistence regions were found for both clusters. Again, the 
coexistence ranges obtained from the molecular dynamics calculations coincide very well 
with the transition ranges obtained from the heat capacities. The lowest energy isomers for 
these clusters represent the building up from A^ = 7 by the successive addition of atoms 
to the five faces of the pentagonal bipyramid core, with the A^ = 12 structure then being 
one atom short of an icosahedron. This regular progression is reflected in the heat capacity 
curves. The A^ = 9 peak is small and close enough to the dissociation peak that it is barely 
more than a shoulder. The A^ = 10 and A^ = 11 peaks are slightly larger, and remarkably 
alike. The A^ = 12 peak is considerably larger than the others. The transition ranges also 
show quite a regular progression, as is evident in Fig. |^, with the A^ = 9 range being the 
smallest and occurring at the lowest temperature. The A^ = 10, 11 and 12 transition ranges 
are about equal in size, but occur at increasingly higher temperatures. 

B. Icosahedral to Double Icosahedral: 13 < A^ < 19 

The results for the range A^ = 13 to 19 are consistent with the molecular dynamics results 
for the magic number endpoints A^ = 13 and 19, but there are some surprising differences for 
sizes in between. A^ = 13 is the prototypical magic number cluster and has been the subject 
of many investigations. Its minimum energy configuration is the icosahedron, a compact 
structure that is extremely stable. The A^ = 13 potential energy surface is characterized 
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by very deep wells, corresponding to the icosahedral isomers, which are separated by large 
barriers and lie far below the wells corresponding to the next lowest energy isomers. Its heat 
capacity behavior is a reflection of the shape of this potential energy surface. The curve 
has a very large, pronounced peak in the solid-liquid region (more than double the height 
of the preceding = 12 peak), and the fluctuations in the heat capacities obtained from 
Metropolis Monte Carlo simulations show the characteristic features of a coexistence region. 
The coexistence range obtained from molecular dynamics calculations reported in Ref. |^ 
also coincides very well with the transition region found in this study. 

The J-walking results also agree with the original J-walking results for Aris reported 
in Ref. except for a small discrepancy in the peak height where a value of {Cy)/N = 
9.37 ± 0.04 was obtained at T* = 0.285 ± 0.004, compared with {Clr)/N = 9.06 ± 0.01 
obtained at Tp = 0.2875 ± 0.0007 in this study. This investigation was more thorough than 
the original one. In that study, computer constraints limited the J-walker distribution sizes 
to only 5 x 10^ configurations, and walk lengths for each temperature consisted of 10'' moves. 
In this study, J-walker distributions were twenty times larger, and walks consisted of 10^ 
passes, and so were thirteen times longer. In addition, the first J-walker distribution in the 
original investigation was obtained from a long Metropolis Monte Carlo walk at T* = 0.42. 
This was a temperature thought to be in the liquid region (and thus presumed to be free 
of quasi-ergodicity), but was actually near the beginning of the dissociation region where 
quasi-ergodicity again begins to become problematic. Since subsequent lower temperature J- 
walker distributions were generated in sequence from this distribution, any systematic errors 
in the distribution could have propagated to the other distributions and thus affected the J- 
walking results. In this study, the initial J-walker distribution was generated at T* = 0.618, 
a temperature on the high temperature side of the dissociation region heat capacity peak 
where the heat capacity fluctuations had dropped substantially, thus reducing systematic 
errors in the lower temperature distributions. Despite these problems, the earlier J-walking 
study still gave a value for the heat capacity peak that was within 4% of the value obtained 
in this study for not much more computational overhead than the similar Metropolis Monte 
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Carlo simulation also reported in Ref. which had an error of 19% for the peak height. 

Further support for the accuracy of the results reported here can be been found in a 
recent study by Tsai and Jordan,ii who combined histogram and J-walking methods to 
calculate heat capacity curves for Aris. They obtained a peak height of {Cy)/N = 9.01 at 
Tp = 0.286, which is very close to the value reported here; the size of the difference between 
my values and theirs is similar to the small variations they observed for results obtained using 
different constraining sphere radii (their results were obtained using a smaller constraining 
radius, with Rc = 10 A, or 2.94a, compared to Rc = 4a used in this study). Finally, the 
results obtained from the Metropolis runs of the same length (10^ passes) agree with the 
J-walker resuhs, with (C^)/A^ = 8.97 ± 0.04 and = 0.2865 ± 0.0011. This indicates that 
Metropolis walk lengths of 10^ passes are sufficiently long to overcome quasi-ergodicity for 

= 13, consistent with the findings of Tsai and Jordan, who reported that walk lengths of 
6 X 10® passes were sufficiently long. 

The results obtained for A^ = 14 to 18 are surprising, and are considerably different 
from the results obtained for A^ = 8 to 12. This was unexpected since the two ranges are 
analogous in many respects. Both are bounded by magic number sizes (A^ = 7 and 13, and 
A^ = 13 and 19, respectively) and have a similar "growth" sequence in their minimum energy 
configurations. This can be seen in Fig. H, which shows the minimum energy configurations 
for A^ = 8 to 19.0 Where the N = 8 lowest energy isomer consists of a lone atom located on 
one of the ten faces of the N = 7 pentagonal bipyramid core, and where the succeeding lowest 
energy isomers can be formally obtained by the addition of atoms in a fivefold symmetric 
manner, with the sixth capping the ring to form the A^ = 13 icosahedron, the A^ = 14 
lowest energy isomer likewise consists of a lone atom located on one of the twenty equivalent 
icosahedral faces, and the succeeding lowest energy isomers (except A^ = 17) can be likewise 
formally obtained by the addition of atoms in a fivefold symmetric manner, with the sixth 
capping the ring to form the A^ = 19 double icosahedron (the expected A^ = 17 configuration 
obtained in this manner is only slightly higher in energy than the actual minimum energy 
isomer shown, with an energy V* = —61.3071, compared to V^j^ = —61.3180). The striking 
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similarity in the binding energy differences for these two ranges can be seen in the plot 
of AEf,{N) shown in Fig. |^ — the sequence for = 8 to 13 is nearly identical to the 
sequence for = 14 to 19. Molecular dynamics calculations reported in Ref. |^ also found 
analogous behavior. Both Atj and Ari3 had relatively large coexistence ranges, while neither 
Arg nor Ari4 showed any coexistence behavior, and both Arg and Aris showed coexistence 
behavior, but had narrower ranges than Ary and Ari3. Wales and Berryll ran extensive 
molecular dynamics simulations combined with quench studies for Aij, Arg, Aria and Ari4 
that confirmed the earlier results. Yet, the heat capacity curves obtained with J-walking 
depicted in Figs. ^ and |^ are not at all consistent with these results. Where the A^ = 8 
heat capacity curve is featureless in the solid-transition region, the A^ = 14 curve shows a 
clear peak. The peak is much diminished relative to the A^ = 13 peak, but is nonetheless 
quite substantial; it is larger than the A^ = 12 peak, which is the largest peak for A^ < 13. 
The fluctuations in the Metropolis heat capacities are consistent with the J-walking results, 
showing only a moderately steep rise followed by a slow decrease, suggesting some coexistence 
behavior. The results for A^ = 15 are even more surprising in comparison. Ref. |^ reported a 
coexistence region for Aris, analogous to Arg, but the J-walking results reveal only a shoulder 
in the heat capacity curve over this region. While this does not rule out a coexistence 
region (as evidenced by the shoulder in the heat capacity curve for N = 7, which does 
show coexistence behavior), it does, then, make the absence of the coexistence region in 
the A^ = 14 molecular dynamics simulations look even more puzzling by comparison. The 
behavior of the Metropolis heat capacity fluctuations for A^ = 15 is also very different; 
the typical sharp rise signifying the onset of the coexistence region has been replaced by 
a more gradual increase that rises slowly over the domain spanned by the shoulder before 
beginning the usual rise in the dissociation region. The inflection points characterizing the 
heat capacity curve occur well after the molecular dynamics coexistence region, and the rise 
in the fluctuations occurs well before, resulting then in a very much larger transition region. 

Rms bond length fluctuation curves for A^ = 13, 14 and 15 are shown in Fig. ^. As with 
N = 7, magic number behavior is manifested in A^ = 13 by its much higher temperature T/ 
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signifying the sharp rise in S{T). Further similarities between = 14 and = 15 are evident 
in the S{T) curves. Except for the unusual smaller rise in 6{T) beginning at T* = 0.064 
for = 14, the two curves differ very little. This smaller rise is likely due to the facile 
movement of the extra atom for = 14 on the icosahedral core. The average potential 
energy corresponding to this temperature is V* = —46.60, which is comparable to the 
saddle point energy of —47.06 representing the edge-bridging isomerization that moves the 
fourteenth atom from one icosahedral face to a neighboring one. It is also comparable with 
two other saddle point energies of —45.95 and —45.84, which also represent isomerizations 
to degenerate capped icosahedrons.i 

A^ = 16 continues the trend from A^ = 15, with the heat capacity curve reduced to 
an even smaller shoulder, and with the Metropolis heat capacity fluctuations remarkably 
uniform throughout the solid-liquid region. The results for A^ = 17, 18 and 19 mirror the 
results for A^ = 13, 14 and 15, with A^ = 17 having a shoulder in its heat capacity curve 
similar to the A^ = 15 curve, A^ = 18 having a moderate peak similar to the A^ = 14 peak, 
and A^ = 19 having pronounced peaks in both the heat capacity and the Metropolis heat 
capacity fluctuation curves, consistent with its magic number status. The only significant 
difference is that the shoulder and peak sizes for A^ = 17, 18 and 19 are smaller than their 
counterparts for A^ = 15, 14 and 13. Interestingly, the rms bond length fluctuation curve 
for A^ = 18 also mirrors the A^ = 14 curve with a smaller rise beginning at T* = 0.106 and 
ending at T* = 0.168, where it rises sharply. Further evidence of magic number behavior in 
A^ = 19 is evident in the S{T) curves shown in Fig. ^ where once again can be seen to 
occur at a much higher temperature than that of the next two larger clusters. The molecular 
dynamics coexistence range for A^ = 19 reported in Ref. |^ also coincides quite well with the 
transition range I obtained, although it is slightly contracted. 
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C. Beyond Double Icosahedral: 20 < iV < 24 

The results for = 20 to 24 also differ substantially from those of = 14 to 18. 
This is not unexpected since the minimum energy configurations for these clusters show a 
considerably different structural sequence, with additional atoms occupying sites between 
the three equatorial 5-atom rings of the = 19 double icosahedron. Fig. [TO shows the 



minimum energy configurations for = 20 to 24. For A^ = 23, the additional atom caps an 
interpenetrating icosahedral substructure, resulting in a magic number configuration. Magic 
number behavior can be seen in Figs. ^ and ^ by both the large heat capacity peak and in 
the binding energy difference maximum. The heat capacity peak is about the same size as 
the previous magic number cluster, A^ = 19. The Metropolis heat capacity fluctuations also 
show the sharp rise and larger peak indicative of a coexistence region. 

Unlike A^ = 14, whose heat capacity peak is much smaller than its preceding magic 
number A^ = 13 peak, the peaks for A^ = 20 and 24 are only slightly smaller than their 
preceding magic number peaks. The major differences between A^ = 19 and 20, as well 
as between A^ = 23 and 24, lie at the beginning of the transition regions. While the heat 
capacity fluctuations for A^ = 19 and A^ = 23 show the sharp rise and subsequent slow 
decline characteristic of a coexistence region, the fluctuations for A^ = 20 and A^ = 24 show 
a more gradual rise to a smaller maximum, followed by a slow decline. This was also the case 
for A^ = 21 and 22, and so determining the onset of the transition region for each of these 
clusters was ambiguous. As with A^ = 14, there was no evidence from molecular dynamics 
calculations of a coexistence region reported for A^ = 20 in Ref. ^, which is surprising since its 
heat capacity peak is nearly as large as that for A^ = 19 and its heat capacity fluctuations are 
also quite pronounced. The heat capacity curves and Metropolis heat capacity fluctuations 
are remarkably alike for A^ = 21 and 22, reminiscent of the similarity between A^ = 10 and 
11. Both A^ = 21 and 22 also have smaller peaks than the other clusters in this range. 
Molecular dynamics results for N = 22 were also reported in Ref. ^ and no evidence of a 
coexistence region was found. 
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IV. DISCUSSION 



Although the heat capacity-temperature curves do not provide the definitive measure 
of the cluster solid-liquid transition region, they do serve as yet another useful probe into 
this region, providing additional data of surprising diversity and complexity that needs to 
be examined in terms of the various theories that have been proposed. Many of the results 
reported here are consistent with the view that for some clusters (especially for the magic 
number sizes), the solid-liquid transition region corresponds to a coexistence region with 
sharply defined, but separate, melting and freezing temperatures. Just how sharp these 
transition temperatures really are remains an open question, but for some clusters, prop- 
erties obtained from Metropolis simulations such as the heat capacity fluctuations and rms 
bond length fluctuations do exhibit threshold-temperature behavior, with abrupt changes 
occurring over narrow temperature ranges. For many of the clusters where molecular dy- 
namics simulations found evidence of coexistence behavior, the Metropolis heat capacity 
fluctuations were consistent, and the transition ranges obtained from Metropolis heat ca- 
pacities coincide quite well with the coexistence ranges. But there are also some puzzling 
discrepancies, such as with = 14 and 15. Molecular dynamics simulations found evi- 
dence of coexistence behavior for N — 15, which has neither a heat capacity peak nor an 
appreciable change in the heat capacity fluctuations throughout the solid-liquid transition 
region, but found no coexistence behavior for N = 14, which does have a heat capacity 
peak and an increase in the fluctuations. On the other hand, the view of the solid-liquid 
transition being a series of isomerizations through an ever increasing number of structures 
also appears consistent with my findings for some of the clusters studied, but it seems hard 
pressed to account for the entire range of behavior seen in the heat capacity curves and 
rms fiuctuations. In addition, some supporting evidence is contrary to the heat capacity 
behavior. For example, the instantaneous normal mode analysis of = 7 and 8 indicated 
similar behavior in the transition region, but the heat capacity studies showed them behav- 
ing very differently. Similarly, the instantaneous normal mode analysis oi N — 13 indicated 
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no abrupt changes for some dynamical properties in the transition region, where the heat 
capacity, its fluctuations, and the rms bond length fluctuations all show abrupt changes. 

It is also clear (for these smaller clusters, at least) that heat capacity peaks by themselves 
are incomplete measures of the solid-liquid transition region. For example, both = 19 
and = 20 have similarly sized heat capacity peaks, but = 19 also shows magic number 
behavior in its binding energy and rms bond length fluctuations, and shows coexistence 
behavior according to molecular dynamics simulations, while = 20 exhibits none of these. 
Even the absence of a heat capacity peak in the solid-liquid region is no indication of the 
absence of a transition there. The curves for A^ = 4, 5 and 8 are all featureless in this 
region, but rms bond length fluctuation curves for A^ = 5 and 8 show evidence of a solid- 
liquid transition. A^ = 7 is a magic number cluster, but its heat capacity curve has no peak, 
even though its non-magic number predecessor A^ = 6 does have a heat capacity peak. 

Additional insight into the solid-liquid transition region can be obtained by comparing 
the heat capacity magic number behavior with that of other properties. Fig. ^ shows that 
the heat capacity magic numbers are the same as the AEi,{N) magic numbers over the 
range 4 < A^ < 24, but there are interesting differences between the two sequences. The 
magic numbers for the AEi,{N) sequence are N = 7, 13, 19 and 23, as evidenced by these 
sizes having the largest values (except for N = 7, which has the largest value for A^ < 9). 
In each case, AEfy{N) substantially decreases for A^ immediately following, and then rises. 
The magic numbers for the heat capacity are likewise A^ = 7, 13, 19 and 23, which are the 
local maxima in the heat capacity peak sequence. However, with the exception of A^ = 8, 
which has a featureless transition region, the peaks for A^ = 14, 20 and 24 do not show a 
corresponding dramatic decline followed by an increase. Instead, there is a gradual decline 
over several sizes. A^ = 14, which has a small peak, is followed by three clusters having 
smaller shoulders, while A^ = 20 has nearly as large a peak as the A^ = 19 peak, and A^ = 24 
has a peak only slightly smaller than the A^ = 23 peak. Furthermore, the heat capacity 
peak sequence for A^ = 8 to 13 does not show the remarkable similarity with the A^ = 14 
to 19 sequence that is evident in the AEb{N) values. Thus, it appears a cluster's minimum 
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energy structure has a dominant influence on the heat capacity only for the magic number 
sizes. 

The impact of the minimum energy structure on cluster energetics for magic number 



sizes is also evident in Fig. [Tl|, which plots for all the clusters studied the potential energies 
per atom as functions of temperature. The curves have a roughly uniform spacing for 
temperatures corresponding to the liquid region, while for temperatures in the solid region, 
the curves show the nonuniform spacing expected for the very different geometrical structures 
associated with each minimum energy configuration, with the enhanced stability of the magic 
number clusters = 13, 19 and 23 (and to a lesser extent, N = 7) clearly evident. Also 
included in the plot are the T* values and the heat capacity peak temperatures, T£, Tp and 
T^, which characterize the transition region. The nonuniform spacing of the solid region and 
the enhanced stability of the magic number structures persist through much of the transition 
region, as does the uniform spacing of the liquid region. This indicates that both solid-like 
and liquid-like forms are indeed found in the region spanned by the heat capacity peak. 
The persistence of solid-like forms into the transition region and the enhanced structural 
stability of the magic number sizes can be seen even more dramatically in Fig. where the 



AEii{N) = —AV*{N) curve for zero temperature shown in Fig. ^ has been generalized to 
nonzero temperatures to form the surface -A{V*{N, T)) = -[{V*{N, T)) - {V*{N-1, T))], 
with the temperature range spanning the solid-liquid transition region up to the beginning 
of the dissociation region of the larger clusters. The icosahedral based magic number sizes 
N = 13, 19 and 23 are especially noticeable — the —A{V*) peak values differ very little 
from their zero temperature values well into the transition region, until their abatement 
over a relatively narrow temperature range at the heat capacity peak temperature. This 



temperature range is much smaller than the transition range defined in Eq. [TT], beginning 
at about T£ rather than at T* and ending well before T^. Although this seems to indicate 
more the abrupt change from solid to liquid rather than the coexistence of solid and liquid 
forms over an extended range, it must be remembered that these energy differences are 
between average energies, which can be quite insensitive to the makeup of their underlying 
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distributions. 

The average potential energy and configurational heat capacity at a given temperature 
are of course just the first two moments of the potential energy distribution for that tem- 
perature, and so they represent only a minimal description of the distribution. Further 
insight into cluster energetics can be obtained by examining the distributions themselves. 
With J-walking, accurate potential energy distributions can be obtained during a simula- 
tion at a given temperature by binning the potential energies during the walk and forming 
a histogram. For example, the distributions shown in Fig. § were obtained this way. By 
projecting several histograms over a temperature range, a distribution surface can be formed 



as a function of the potential energy and temperature. Fig. |1^ shows examples of such sur- 
faces. These are for the range = 13 to 18, which contains representative samples of the 
different types of heat capacity behavior encountered (a strong peak at = 13 declining 
to a weak shoulder at = 16 and then rising again to a strong peak at = 19). Narrow 
distributions occur at the temperature extremes corresponding to the solid and dissociation 
regions; saddle points correspond to the widest distributions and thus represent heat ca- 
pacity peaks. Not surprisingly, the change in the distributions in the sohd-liquid transition 
region is much more abrupt for A^ = 13 than for the other clusters. Not only do the distribu- 
tions for the solid region widen appreciably and then narrow suddenly for the liquid region 
to form the characteristic hump on the surface, but the distributions for the two regions are 
partially offset, indicating a bimodal distribution. This is evident in Fig. |14|, which shows 
the potential energy distributions for temperatures near the heat capacity peak temperature 
Tp. The bimodal distribution indicates well separated energies that can be associated with 
distinct solid-like and liquid-like forms and thus are further evidence of coexistence behavior 
in A^ = 13. Unfortunately, A^ = 13 is the only cluster in the range studied that clearly has a 
bimodal distribution in the solid-liquid transition region. As can be seen in Fig. 0, even the 
other magic number clusters A^ = 19 and A^ = 23, which have similar distribution surfaces 
with a similar widening at the heat capacity peak temperature, do not show clear bimodal- 
ity there. This does not rule out coexistence behavior, though, since the two distributions 
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might overlap enough to show only a single maximum, but it does indicate that the different 
distributions are not widely spaced apart, like they are in = 13. 

The distribution surfaces for = 15, 16 and 17 show a very smooth, gradual evolution 
from the sharp, narrow distributions in the solid region to the broader distributions in 
the liquid region. This behavior is consistent with the view of the solid-liquid transition 
being a progression through an increasingly large number of structural isomers. For these 
clusters then, there are no threshold changes in the potential energies being accessed as the 
temperature increases past the solid region — higher energy isomers with energies slightly 
higher than the minimum energy isomer are separated by low barriers and thus are easily 
accessible. = 16 clusters have several local minima lying close to the global minimum,0 as 
does = 17,0 which also has a local minimum differing by less than 0.02%. The similarity 
between these three clusters is also evident in the behavior of their average potential energy 



per atom. As can be seen in Fig. |TT|, the spacings between the curves for these clusters are 
almost equal and change very little throughout the temperature range. Again, this is very 
much different than the behavior displayed by the magic number clusters. 

Further evidence of the coexistence between solid-like and liquid-like isomers for the 
magic number clusters in the transition region comes from the rms bond length fluctuation 
curves. Comparison of the S(T) curves shown in Fig. ^ with the —A{V*{N,T)) surface in 
Fig. reveals that the temperatures where the fluctuations for each cluster rise sharply 
from their nearly linear increase in the solid region (thus marking the onset of hindered 
isomerization) are well below the temperatures where the —A{V*{N,T)) peaks diminish 
quickly, implying that the solid-like forms persist well into the transition region (although 
undergoing increasingly frequent isomerizations) . However, the end of the sharp rise and 
the leveling of the S{T) curves to a value of about 0.3 also occurs in this region, indicating 
liquid-like behavior there as well. Note also that the lack of change in the average potential 
energy difference between A^ = 14 and A^ = 13 up to a reduced temperature of about 0.25 
lends further support to the conjecture that the initial rise in S{T) for A^ = 14 is due to the 
motion of the lone atom on the icosahedral core, rather than an isomerization to one of the 
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other low-lying configurations. 



Magic number behavior in Tg is shown in Fig. 15, which plots the values obtained 



from the Metropolis simulations against cluster size. Magic number behavior is evident for 
= 13, 19 and 23 (and to a lesser extent for = 7) by their much higher values. 
The values for all the clusters lie below T* = 0.18, which according to Figs. |TI| and ^ 
places them near the end of the solid region. Also, as can be seen in Fig. |l^, they generally 
coincide with the T* values obtained from the rise in the heat capacity fluctuations, further 
supporting the claim that this region marks the beginning of the transition. Given the large 
degree of ambiguity in determining either T* or for some of the larger clusters because of 
the large levels of noise in their fluctuations, it is not valid to conclude much more than the 
general observation that these values increase slightly with increasing cluster size, except for 
magic number sizes, where they occur at much higher temperatures. 

Also included in Fig. |15| are the T^j^^j^ values obtained from molecular dynamics sim- 
ulations reported in Ref. |^. Although these values are all significantly higher than the 
Metropolis values, their pattern is quite similar, with the magic number clusters having val- 
ues that are substantially higher than those of their immediate neighbors. As noted earlier, 
the discrepancies between the two sets of data result from the shorter trajectory lengths 
used in the molecular dynamics simulations, and from the different ensembles associated 
with each method.S 

Those clusters exhibiting coexistence behavior also have their coexistence ranges shown 
in Fig. |15|. Except for A^ = 15, which has a coexistence range, contrary to the heat capacity 
behavior, and A^ = 20, which has no coexistence range even though the heat capacity data 
suggests it should, the correspondence between the coexistence ranges and the transition 
ranges defined by the heat capacity behavior is quite good. For A^ = 7, 9, 11, 13 and 19, 
both studies give consistent evidence of coexistence behavior, and the size and location of 
the two ranges are similar, while for A^ = 8, 17 and 22, the lack of evidence of coexis- 
tence behavior is likewise consistent. Unlike most of the molecular dynamics results used 
as evidence for the coexistence region, though, the heat capacity results obtained from the 
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Monte Carlo simulations are equilibrium properties and thus provide mostly inferential, not 
direct, evidence of dynamical coexistence between solid-like and liquid-like forms. The large 
degree of consistency between the two does lend strong support to the hypothesis, but the 
discrepancies for = 15 and = 20 are troubling in that they have no obvious resolu- 
tion. While it appears that quasi-ergodicity is an even more insidious problem than has 
been appreciated in the past, one cannot simply dismiss these discrepancies as being due to 
improper sampling in the molecular dynamics simulations without accounting for the good 
agreement for the other sizes, which were simulated in similar manner (also, it is the magic 
number sizes that are most prone to quasi-ergodicity problems, and these are the sizes where 
agreement is the strongest). The discrepancies might be the result of the different ensem- 
bles used,@ since some cluster properties have a strong ensemble dependence. However, 
comparison of ensemble sensitive properties such as the rms bond length fluctuations and 
diffusion constants that have been calculated in both the micro-canonical and canonical en- 
sembles reveals more similarities than differences (the the temperature curves have roughly 
the same shape and differ primarily by their displacement in temperature) and again, one is 
left with the dilemma of justifying why the agreement between the Monte Carlo results and 
the molecular dyanamics results is so poor for these clusters but so good for the other sizes. 

Despite these inconsistencies, the bulk of my results are consistent with the view of a 
dynamic coexistence in the transition region and can be interpreted within that framework. 
My results are not totally inconsistent with the viewpoint that emphasizes the smooth 
progression of isomerizations between locally similar but globally distinct configurations 
occurring throughout the transition region. Clearly, cluster transitions are a progression of 
isomerizations — the issue is just how smooth these progressions really are. For clusters such 
as = 15, 16 and 17 the progressions appear to be very smooth, but for many other clusters 
(especially the magic number clusters), the threshold behavior of many of their properties 
and their very different heat capacity behavior implies an abruptness in the progressions 
that is more appropriately described then by the coexistence viewpoint. 

The nature of cluster phase transitions still remains much unresolved then, despite the 
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additional information provided by the heat capacity results. More data is required to resolve 
the outstanding issues. This can be accomplished by extending the systematic survey of 
cluster properties beyond N = 24, but another avenue also seems appropriate. Since many 
cluster properties result from an interplay of particle dynamics based on coordinated atom 
moves, and structural effects based on particle size and intermolecular forces, the study 
of heterogeneous atomic clusters of various size and composition can help determine the 
relative importance of each factor. Work is currently underway on a J-walking study of 
binary Ne-Ar and Ar-Kr clusters. Ne and Ar differ greatly in both size and potential, but 
Ar and Kr have similar sizes and so differences in their behavior can provide more insight 
into the transition region. 
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23 were performed on a DEC Alpha 3000/300L AXP equipped with 32 Mb of RAM. 

The Savitsky-Golay algorithm was used for the data smoothing, and 2D cubic splines were 
used for the interpolation. These routines are implemented in the MS-DOS graphics soft- 
ware package Axum, version 3.0. Different fits were obtained by using different smoothing 
orders, corresponding to 7, 9, 11 or 13 surrounding data points being used to fit each 
point. 

The self diffusion constants as functions of temperature differed substantially from those 
obtained from molecular dynamics calculations by Beck and Marchioro (Ref. ^ by their 
lack of a sharp rise in the transition region. A subsequent instantaneous normal mode 
analysis by Adams and Stratt (Ref. ^ using an expression for the diffusion constant that 
was more sensitive to the low frequency part of the phonon spectrum gave results that 
were in agreement with the molecular dynamics calculations. The curves for Atj and Arg 
remained very much alike. 

Because the J- walking runs were started at high temperatures from random configurations, 
and the temperatures were decreased in small increments down to zero, the simulations 
also served as global minimization calculations, akin to simulated annealing. In all cases, 
the lowest energy configuration obtained with J-walking was the same as that reported in 



Refs. 29 and 31. 
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FIGURES 

FIG. 1. Classical internal energy (at left) and constant volume heat capacity (at right) as func- 
tions of temperature for Ar23 clusters; both are in reduced units, with U* = U/e and Cy = Cv/ks- 
The lower temperature heat capacity peak corresponds to a solid-liquid transition, while the higher 
temperature peak corresponds to a liquid-vapor transition. The solid curves were obtained using 
J-walking from external J-walker distributions containing 10^ configurations; the jump attempt fre- 
quency was Pj = 0.1. For each temperature, 10^ warmup passes were run, followed by 10'' passes 
with data accumulation. The dashed curves were obtained using standard Metropolis Monte Carlo 
methods, likewise consisting of 10^ warmup passes followed by lO'' passes with data accumulation. 
Some representative single standard deviation error bars have been included. Also included in the 
heat capacity plot are results from Metropolis runs consisting of 10^ warmup passes followed by 10^ 
passes with data accumulation (curve i) and 10^ passes with data accumulation (curve ii). Lack of 
convergence in the Metropolis simulations is evident in the solid-liquid transition region. 

FIG. 2. Normalized potential energy histograms for = 11 and A = 19 J-walker distributions. 
Distribution particulars are given in Tables | and As the cluster size increases, so does the 
potential energy range, implying more distributions are required to span the entire range. 

FIG. 3. The dependence on walk length for Metropolis averages of the rms bond length fluc- 
tuations 5 as a function of temperature for A^ = 13. The values at each temperature are averages 
obtained from 10 Metropolis Monte Carlo walks having lengths of lO'' passes (short dashes), 10^ 
passes (long dashes) and 10^ passes (solid line). The dotted curve is the heat capacity from Fig. ^, 
and has been included for comparison. Slow convergence due to quasi-ergodicity in the transition 
region is evident by the displacement of the curves to lower temperatures with increased walk 
length. The insert shows the temperature = 0.173 where the rms bond length fluctuations be- 
gin to rise sharply. The solid triangle indicates the temperature T^md where the rms bond length 
fluctuations obtained from molecular dynamics simulations reported in Ref. ^ were found to rise 
sharply. 
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FIG. 4. Heat capacity curves as functions of temperature (upper curves in each plot, in re- 
duced units) for clusters ranging from = 4 to 12. All curves where obtained using J- walking 
from externally stored J-walker distributions with jumps attempted with a frequency of Pj = 0.1. 
The distribution parameters are listed in Table |. For each temperature, simulations consisted of 
10^ warmup passes followed by lO"^ passes with data accumulation. The lower, noisy curves are 
heat capacity standard deviations obtained from similar Metropolis Monte Carlo simulations, each 
consisting of 10^ warmup passes followed by 10'' passes with data accumulation. The tempera- 
tures corresponding to a sudden increase in standard deviation (indicated by the arrows) serve as 
estimates to the beginning of the solid-liquid transition regions. For those clusters having heat 
capacity curves with peaks in the solid-liquid transition region (A^ = 6, 9-12), the temperatures 
corresponding to the low temperature side inflection point, T£, the peak, Tp, and the high temper- 
ature side local minimum, are indicated on the temperature axis by the large ticks. For N = 7, 
which has a shoulder in its heat capacity curve in this region, the large ticks indicate the inflection 
points, T£ and T^. These peak parameters are also listed in Table For those clusters having 
molecular dynamics simulations reported in Ref. |3| {N = 7-9 and 11), the solid triangles indicate 
the temperatures at which the rms bond length fluctuations 6 were found to rise sharply, while 
the horizontal error bars represent the coexistence temperature ranges, AT* = — TJ, for those 
clusters found to have a bimodal form in the distribution of short-time averaged kinetic energies. 

FIG. 5. Same as Fig. |^, but for clusters ranging in size from = 13 to 24. The J-walker 
distribution parameters are listed in Table ||. 
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FIG. 6. Magic number effects in cluster heat capacities. The top plot shows the solid-liquid 
transition region heat capacity peak characteristic temperatures as functions of cluster size A^. For 
those clusters having distinct peaks, the values T£, Tp and are indicated, while for those clusters 
having heat capacity shoulders in this region, Tp is taken to be T^. These points are also indicated 
in Figs. ^ and ^ by the large ticks on the temperature axes. No solid-liquid transition region heat 
capacity peaks or shoulders were found for = 4, 5 and 8. The lower T* curve indicates the 
dependence on cluster size for the beginning of the transition region (also indicated by the arrows 
in Figs. ^ and ^). The solid-liquid transition ranges can be defined as the difference between the 
top and the bottom curves, AT* = — T*. The middle plot shows the size dependence on the 



corresponding heat capacities. The bottom plot was inspired by Ref. 29 and has been included 
for comparison. It plots the binding energy differences between successive clusters, as a function 
of size, AEb{N) = —AVmm{]^) = — [V'minl-^) " Vinm{]^ ~ !)]• Magic number sizes based on zero 
temperature structural stability (N = 13, 19 and 23, and to a lesser extent, N = 7) have a clear 
correlation with those evident in the heat capacity peaks. 

FIG. 7. Rms bond length fluctuations as functions of temperature for N = 4-9. The values at 
each temperature are averages obtained from 10 Metropolis Monte Carlo walks having lengths of 
10^ passes. Some representative single standard deviation error bars have been included. Only the 

= 4 curve shows no abrupt change in the solid-liquid region before rising sharply in the dissoci- 
ation region. The N = 7 curve shows magic number behavior that reflects the enhanced stability 
of the cluster's lowest energy conflguration. Its sharp rise occurs at a much higher temperature 
than that for either = 6 or = 8. 
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FIG. 8. Pentagonal "auf-bau" in the minimum energy configurations for Lennard-Jones clus- 
ters. The configurations shown are the lowest energy isomers obtained from the J-walking sim- 
ulations, and are identical to those found in other studies.ii~0 Beginning with the pentagonal 
bipyramid core corresponding to = 7, the minimum energy configuration for each succeeding 
cluster up to = 12 can be formally obtained by adding an atom (indicated by the shading) to 
a neighboring pentagonal site to form a second five-atom ring, which is then capped to form the 

= 13 icosahedron. Except for N = 17, the minimum energy configurations for = 14 to 19 
are analogous, with succeeding clusters built up from the icosahedral core by adding an atom to 
a neighboring pentagonal site to form a third five-atom ring, which is then capped to form the 
A^ = 19 double icosahedron. This correspondence between clusters A^ = 8 to 13 and A^ = 14 to 19 
is reflected in the binding energy differences, AEj,, shown in Fig. ^ — the sequence for A^ = 8 to 
13 is nearly identical with the sequence for A^ = 14 to 19. 

FIG. 9. Rms bond length fluctuations as functions of temperature for A^ = 13-15 and 
A^ = 19-21. The values were obtained in the same manner as those in Fig. |^. Magic number 
behavior is also evident for A^ = 13 and A^ = 19 by their significantly higher temperatures marking 
the onset of the sharp rise in S(T). The A^ = 14 curve is somewhat unusual, showing a smaller, 
nearly linear rise at about T* = 0.064 until about T* = 0.13 where it then undergoes the sharp 
rise typical of the other clusters. 

FIG. 10. Minimum energy configurations for A^ = 20 to 24. The configurations shown are the 
lowest energy isomers obtained from the J-walking simulations, and are identical to those found 
in other studies.@'0 Each configuration can be obtained formally from the preceding one by the 
addition of the shaded atom. For A^ = 23, the additional atom caps an interpenetrating icosahedral 
substructure, forming a magic number structure. 
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FIG. 11. Potential energy curves as functions of temperature obtained from J- walking simula- 
tions for clusters ranging from = 4 (top curve) to = 24 (bottom curve). The magic number 
sizes {N = 7, 13, 19 and 23) are indicated by the dashed curves. Reduced units have been used, 
with V* = V/e and T* = ksT/e. The dotted lines intersecting the curves are the values for the 
beginning of the transition region T* (as characterized by a sharp rise in the heat capacity fluctu- 
ations), and the solid- liquid transition region heat capacity peak characteristic temperatures (T£, 
Tp and for clusters having distinct peaks, or the inflection points T£ and for those clusters 
having heat capacity shoulders). The transformation from the solid-like region to the liquid- like 
region is evident by the change in spacing between the curves from nonuniform at low temperatures 
to nearly uniform at high temperatures. 

FIG. 12. Potential energy differences as functions of temperature T* and size (the vertical 
axis is -A(y(A^,T)) = -[(VAr(r)) - {Vn-i{T))]). The temperature mesh size is AT = 0.2 K, 
or AT* = 0.00562. The structure of the zero temperature curve AEf,(N) (which is also shown 
in Fig. ^) can be seen to extend well into the transition regions, indicating the persistence of 
the solid-like lowest energy isomer, especially for the magic number sizes. The dashed curves 
correspond to the temperatures closest to the heat capacity peak temperatures Tp for the magic 
numbers = 13 (long dashes), 19 (medium dashes) and 23 (short dashes). 

FIG. 13. Normalized potential energy distribution surfaces as functions of potential energy and 
temperature for A^ = 13 to 18. The surfaces are comprised of potential energy histograms obtained 
from J-walking simulations with temperature increments of AT = 0.2 K (AT* = 0.00562). In 
each case, the truncated peak on the left represents the solid region and the peak on the right the 
dissociation region; saddle points correspond to heat capacity peaks. 
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FIG. 14. Normalized potential energy histograms for temperatures near the sohd-liquid transi- 
tion region heat capacity peak for the magic number clusters. The solid curves are for temperatures 
below Tp, the dotted curves for temperatures above; the temperature increment is AT = 0.1 K 
(AT* = 0.00281). The data was obtained from the J-walking simulations. Only the = 13 
distributions have a bimodal form indicative of solid-liquid coexistence. 

FIG. 15. Magic number effects in cluster rms bond length fluctuations. The plot shows the 
temperatures T^*mC' where S{T) obtained from Metropolis simulations begins to deviate sharply 
from the nearly linear increase associated with the solid region (see Fig. |^ for an example). For 
N = 14, 18, 22 and 24, the lower triangle signifies the temperature where S{T) first has a moderate 
rise at the end of the solid region, and the upper triangle signifies temperature where S{T) then 
rises sharply. Also shown are the corresponding values Tj md obtained from molecular dynamics 
simulations reported in Ref. ^; the vertical error bars indicate the coexistence ranges found in the 
molecular dynamics studies. The Metropolis values occur at much lower temperatures, but show a 
similar pattern as a function of size. Also included are the heat capacity peak temperatures Tp and 
the transition ranges AT* = T^ — T*, which are indicated by the dotted lines. The temperatures 
^(S*MC seen to roughly correspond to the beginning of the transition ranges T* for most 

sizes. 
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TABLES 

TABLE L J-walker distributions for 4 < < 12. The distributions for each cluster were 
generated in stages, with the initial J-walker distribution obtained from a single Metropolis walk 
at a temperature high enough to be in the cluster dissociation region, and subsequent lower tem- 
perature distributions obtained using J-walking from the preceding distribution; the jump attempt 
frequency was Pj = 0.1 throughout. For each distribution, 10^ warmup moves were made before 
configurations were stored. All distributions consisted of 10^ configurations stored in 10 or 50 files 
(depending on the amount of computer primary memory), except for A^ = 4 and N = 8, which had 
1.25 X 10^ and 0.75 x 10® configurations, respectively. Configurations were stored periodically at 
intervals of Jextra moves. For those distributions obtained using J-walking, %Ja gives the percent- 
age of attempted jumps that were accepted. Absolute temperatures are for Ne. Potential energy 
histograms for the A^ = 11 distributions are shown in Figure 0. 
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0.112 


600 


8.6 




N = 


7 






A^ = 


8 






A^ = 


9 




r(K) 




J 2xtra 


%Ja 


T(K) 




J 2xtra 


%Ja 


T(K) 




J extra 


%Ja 


15.0 


0.421 


1000 




20.0 


0.562 


800 




20.0 


0.562 


900 




12.0 


0.337 


1000 


18.2 


15.0 


0.421 


800 


45.9 


15.0 


0.421 


900 


28.1 


8.0 


0.225 


1000 


9.1 


12.0 


0.337 


800 


9.3 


11.5 


0.323 


900 


7.1 










8.0 


0.225 


800 


11.1 


7.0 


0.197 


900 


6.2 




N = 


10 






A^ = 


11 






A^ = 


12 
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r(K) T* Jextra %Ja T(K) T* Jextra %Ja T(K) T* Jextra %Ja 

20.0 0.562 1000 20.0 0.562 1000 20.0 0.562 1200 

16.0 0.449 1000 34.2 16.0 0.449 1100 19.8 17.0 0.478 1200 30.8 

12.0 0.337 1000 4.4 12.0 0.337 1100 6.2 14.0 0.393 1200 13.1 

8.0 0.225 1000 9.6 8.0 0.225 1100 7.6 10.0 0.281 1200 7.2 

5.0 0.140 1000 6.5 5.0 0.140 1100 5.9 7.0 0.197 1200 8.6 

4.0 0.112 1200 6.5 
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TABLE II. J-walker distributions for 13 < A'^ < 24. The distributions were obtained as de- 
scribed in Table |. Potential energy histograms for the = 19 distributions are shown in Figure 





N = 


13 






N = 


14 






TV = 


15 






TV = 


16 




r(K) 


T* 


■^cxtra 


%Ja 


T(K) 


T* 


J extra 


%Ja 


T(K) 


T* 


extra 


%Ja 


T(K) 


T* 


J extra 


%Ja 


22.0 


0.618 


1300 




22.0 


0.618 


1400 




25.0 


0.702 


1500 




25.0 


0.702 


800 




18.0 


0.506 


1300 


28.2 


18.0 


0.506 


1400 


18.3 


21.0 


0.590 


1500 


49.6 


20.0 


0.562 


1600 


23.0 


14.0 


0.393 


1300 


5.0 


14.0 


0.393 


1400 


6.0 


17.5 


0.492 


1500 


11.6 


16.5 


0.463 


1600 


9.0 


10.5 


0.295 


1300 


8.4 


10.0 


0.281 


1400 


3.8 


14.0 


0.393 


1500 


11.2 


13.0 


0.463 


1600 


11.5 


7.0 


0.197 


1300 


4.3 


6.0 


0.169 


1400 


2.8 


11.0 


0.309 


1500 


13.5 


10.0 


0.281 


800 


12.3 


4.0 


0.112 


130 


6.5 


3.0 


0.084 


140 


2.0 


8.0 


0.225 


150 


8.5 


7.0 


0.197 


400 


6.0 


2.0 


0.056 


130 


3.6 










6.0 


0.169 


150 


20.5 


4.5 


0.126 


160 


7.3 


















4.0 


0.112 


150 


13.5 


2.5 


0.070 


160 


3.9 


















2.0 


0.056 


150 


2.2 












TV = 


17 






TV = 


18 






N = 


19 






TV = 


20 




T(K) 


T* 


"^extra 


%Ja 


T(K) 


T* 


J extra 


%Ja 


T{K) 




extra 


%Ja 


T(K) 


T* 


J extra 


%Ja 


25.0 


0.702 


850 




25.0 


0.702 


900 




25.0 


0.702 


950 




25.0 


0.702 


1000 




21.0 


0.590 


1700 


33.9 


21.0 


0.590 


1800 


26.0 


21.0 


0.590 


1900 


20.8 


21.0 


0.590 


2000 


16.6 


18.0 


0.506 


1700 


15.3 


18.0 


0.506 


1800 


15.8 


17.5 


0.492 


1900 


10.7 


17.5 


0.492 


2000 


11.3 


15.0 


0.421 


1700 


17.7 


15.0 


0.421 


1800 


17.8 


14.0 


0.393 


1900 


11.6 


14.0 


0.393 


2000 


11.3 


12.0 


0.337 


850 


16.4 


12.0 


0.337 


900 


15.8 


11.0 


0.309 


1900 


12.2 


11.0 


0.309 


2000 


11.2 


9.0 


0.253 


850 


9.3 


9.0 


0.253 


900 


7.6 


9.0 


0.253 


950 


14.0 


9.0 


0.253 


1000 


13.9 


6.5 


0.183 


170 


9.5 


6.5 


0.183 


900 


9.2 


6.5 


0.183 


950 


9.3 


6.5 


0.183 


1000 


9.6 


4.5 


0.126 


170 


13.4 


4.0 


0.112 


180 


4.6 


4.0 


0.112 


190 


5.1 


4.5 


0.126 


200 


11.9 


3.0 


0.084 


170 


13.4 


2.5 


0.070 


180 


7.1 


2.5 


0.070 


190 


7.4 


3.0 


0.084 


200 


10.7 


1.5 


0.042 


170 


1.5 


1.5 


0.042 


180 


6.0 


1.5 


0.042 


190 


6.0 


2.0 


0.056 


200 


11.9 


























1.0 


0.028 


200 


1.0 




N = 


21 






N = 


22 






TV = 


23 






TV = 


24 




T(K) 


T* 


■/extra 


%Ja 


T(K) 


T* 


t/extra 


%Ja 


T{K) 


T* 


i/extra 


%Ja 


T(K) 


T* 


J extra 


%Ja 


25.0 


0.702 


1050 




25.0 


0.702 


1100 




25.0 


0.702 


1150 




25.0 


0.702 


1200 




Zl .u 


u. oyu 






91 n 


u.oyu 


990(1 
ZZUU 


1 1 Q 


91 n 
z± .u 


u.oyu 


ZoUU 


ins 


99 n 
zz .u 


n fii s 

U.Oio 


9zinn 
Z'yruu 


ZO. / 


17.5 


0.492 


2100 


11.7 


17.5 


0.492 


2200 


11.8 


18.0 


0.506 


2300 


18.3 


19.0 


0.534 


2400 


18.7 


14.0 


0.393 


2100 


10.9 


14.0 


0.393 


2200 


10.5 


15.0 


0.421 


2300 


17.1 


16.0 


0.449 


2400 


16.9 


11.0 


0.309 


2100 


10.0 


11.0 


0.309 


1100 


9.6 


12.5 


0.351 


2300 


20.2 


13.0 


0.365 


2400 


14.0 


8.5 


0.239 


1050 


8.4 


8.5 


0.239 


1100 


7.4 


10.5 


0.295 


2300 


17.2 


10.5 


0.295 


2400 


11.3 


6.0 


0.169 


1050 


7.9 


6.0 


0.169 


1100 


7.4 


8.5 


0.239 


1150 


13.0 


8.5 


0.239 


1200 


12.5 


4.0 


0.112 


210 


8.1 


4.0 


0.112 


220 


7.3 


6.5 


0.183 


230 


16.4 


6.5 


0.183 


1200 


15.0 


2.5 


0.070 


210 


5.9 


2.5 


0.070 


220 


5.1 


4.5 


0.126 


230 


9.9 


4.5 


0.126 


240 


8.7 


1.5 


0.042 


210 


4.6 


1.5 


0.042 


220 


4.1 


3.0 


0.084 


230 


8.5 


3.0 


0.084 


240 


7.7 
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2.0 0.056 230 9.3 2.0 0.056 240 8.6 

1.0 n.02S 230 0.5 1.0 0.02S 210 0.1 
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TABLE III. Solid-liquid transition region heat capacity peak parameters for 4 < < 24. 
These values were obtained by smoothing and interpolating the J-walking data shown in Figs. § 
and p. The points (Tp, {Cy)p/N) are the heat capacity peaks, while the points (T^, {Cy)H/N) 
are the local minima occurring on the high temperature side of each peak. For those clusters 
having heat capacity shoulders instead of peaks {N = 7, 15-17), the points {T^ , (Cy) h / are 
the inflection points. The points (T£, {Cy)L/N) are the inflection points on the low temperature 
side of the peak. The uncertainty estimates are standard deviations of the values obtained using 



different curve smoothing parameters. 



IN 












/ AT 










/ AT 








/ AT 


4 
5 
6 


0. 


0652 


±0. 


.0004 


3.112 


±0.008 





.1036 


±0.0008 


3.595 


±0.001 


0.1344 


±0.0003 


3.540 


±0.002 


7 

8 


0. 


125 


±0. 


.002 


3.61 


±0.04 












0.190 


±0.004 


4.39 


±0.01 


9 


0, 


137 


±0, 


.002 


3.81 


±0.04 


0, 


.1943 


±0.0006 


4.470 


±0.002 


0.208 


±0.001 


4.455 


±0.002 


10 


0, 


145 


±0, 


.002 


3.99 


±0.06 


0, 


.1942 


±0.0008 


4.834 


±0.001 


0.240 


±0.002 


4.687 


±0.001 


11 


0, 


162 


±0, 


.002 


4.18 


±0.05 


0, 


.215 


±0.001 


5.010 


±0.002 


0.255 


±0.001 


4.925 


±0.003 


12 


0, 


,224 


±0, 


.004 


5.0 


±0.2 


0, 


.272 


±0.002 


6.326 


±0.003 


0.317 


±0.001 


6.146 


±0.001 


13 


0, 


,2593 


±0, 


.0006 


7.10 


±0.06 


0, 


.2875 


±0.0007 


9.06 


±0.01 


0.3553 


±0.0008 


6.794 


±0.004 


14 


0, 


,263 


±0, 


.002 


5.8 


±0.1 


0, 


.307 


±0.001 


7.124 


±0.002 


0.356 


±0.001 


6.770 


±0.002 


15 


0, 


,265 


±0, 


.002 


5.52 


±0.04 












0.339 


±0.002 


6.589 


±0.007 


16 


0, 


194 


±0, 


.001 


4.29 


±0.03 












0.321 


±0.003 


5.86 


±0.02 


17 


0, 


,209 


±0, 


.004 


4.5 


±0.1 












0.295 


±0.002 


5.538 


±0.006 


18 


0, 


,224 


±0, 


.002 


4.80 


±0.06 


0, 


.274 


±0.001 


5.834 


±0.001 


0.326 


±0.002 


5.652 


±0.002 


19 


0, 


,243 


±0, 


.001 


5.92 


±0.08 


0, 


.2739 


±0.0006 


7.318 


±0.002 


0.352 


±0.002 


5.794 


±0.005 


20 


0, 


,2519 


±0, 


.0004 


5.77 


±0.02 


0, 


.2858 


±0.0003 


7.085 


±0.001 


0.3610 


±0.0008 


5.756 


±0.001 


21 


0, 


,257 


±0, 


.002 


5.37 


±0.06 


0, 


.300 


±0.002 


6.383 


±0.005 


0.371 


±0.001 


5.843 


±0.001 


22 


0, 


,258 


±0, 


.003 


5.5 


±0.1 


0, 


.298 


±0.001 


6.333 


±0.002 


0.3727 


±0.0009 


5.726 


±0.001 


23 


0, 


,265 


±0, 


.001 


6.05 


±0.07 


0, 


.297 


±0.001 


7.383 


±0.003 


0.386 


±0.001 


5.764 


±0.001 


24 


0. 


,257 


±0, 


.001 


5.56 


±0.06 


0, 


.2917 


±0.0002 


6.703 


±0.001 


0.3787 


±0.0006 


5.557 


±0.007 
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(*A)d 




N 



20 21 22 23 24 



0.15 



0.10 



0.05 



0.00 



0, 



0.10 



0.05 



0.00 



0, 



0.10 



0.05 



0.00 




